function pdistmin = v(xt0, yt0)

  tfinal = 100;
  time_steps = 10000;
  t = [1:time_steps]./time_steps .* tfinal;

  n = 3;
  dim = 2;


  xIn = zeros(n, dim);
  pIn = zeros(n, dim);


  % for ni=1:n
  %   for di = 1:dim
  %     xIn(ni, di) = (-1).^ni;
  %     pIn(ni, di) = (-1).^ni;
  %   end;
  % end

  xIn(1, 1) = -1;
  xIn(2, 1) = 1;
  xIn(3, 1) = 0;


  pIn(1, 1) = xt0;
  pIn(1, 2) = yt0;
  pIn(2, 1) = pIn(1, 1);
  pIn(2, 2) = pIn(1, 2);
  pIn(3, 1) = -2.*pIn(1, 1);
  pIn(3, 2) = -2.*pIn(1, 2);



  uIn = reshape([xIn; pIn], 2*n*dim, 1);

 
  

  starttime = clock;
  timeouttime = 15;

  timeoutflag = 0;
  

  opts = odeset('RelTol', 1e-6, 'AbsTol', 1e-9, 'OutputFcn', @(t,u,flag) timeout());
  [tOut, uOut] = ode45(@(t, u) eom(t,u), t, uIn, opts);
  if size(uOut, 1) ~= time_steps
    if timeoutflag == 1
      pdistmin = 0; % timeout
    else
      pdistmin = -1; % acc
    end
    return

    % starttime = clock;
    % [tOut, uOut] = ode15s(@(t, u) eom(t,u), t, uIn, opts);
    % if size(uOut, 1) ~= time_steps
    %   pdistmin = 5
    %   return
    % end
  end
  uOut = reshape(uOut, time_steps, 2*n, dim);
  x = uOut(:, 1:n, :);
  p = uOut(:, n+1:2*n, :);

  pdist = phasedist(x, p, xIn, pIn);
  for i=1:time_steps
    if pdist(i) > 5
      break
    end
  end
  pdist = pdist(i:time_steps);
  pdistmin = min(pdist);


  return

  color = ['b', 'r', 'g'];
  close all
  figure 
  hold on
  for i=1:n
    plot(reshape(x(:, i, 1), time_steps, 1), reshape(x(:,i,2), time_steps, 1), color(i));
  end

  figure 
  hold off
  plot(pdist)

  function stop = timeout()
    stop = 0;
    if etime(clock, starttime) > timeouttime
      stop = 1;
      timeoutflag = 1;
    end

  end

  function pdist = phasedist(x, p, xin, pin)
    for t =1:time_steps
      x(t,:,:) = reshape(x(t,:,:), n, dim) - xin;
      p(t,:,:) = reshape(p(t,:,:), n, dim) - pin;
    end
    x = x.^2;
    p = p.^2;
    pdist = sum(sum(x, 2),3) + sum(sum(p,2),3);
    
  end

  function ut = eom(t, u)
    u = reshape(u, 2*n, dim);
    x = u(1:n,:);
    p = u(n+1:2*n, :);

    f = zeros(n, n, dim); 
    if etime(clock, starttime) < timeouttime
      for i = 1:n
        for j=i+1:n
          dist = 0;
          xdiff = x(i, :) - x(j, :);
          dist = norm(xdiff);
          f(i,j,:) = xdiff./(dist.^3);
        end
      end
      for i = 1:n
        for j=i+1:n
          f(j,i,:) = -f(i, j, :);
        end
      end
    end

    pt = sum(f, 1);



    ut = zeros(2*n, dim);
    ut(1:n,:) = p;
    ut(n+1:2*n, :) = pt;

    ut = reshape(ut, 2.*n.*dim, 1);
  end




end
